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The Effect of Radiative Feedback on Bondi— Hoyle Flow 
around a Massive Star 



m 
o 
o 

o 

Q 

(N 



> 

(N 
m 
(N 

o 

O 



Richard Edgar^* and Cathie Clarke^ 

^Institute of Astronomy, Madingley Road, Cambridge CBS OHA, Great Britain 



2 February 2008 



ABSTRACT 

We apply an alKorit hm for radiative feedback on a dusty flow (detailed in 
lEdear fc Clar ke" ('2003')) to the problem of Bondi-Hoyle accretion. This calculation 
is potentially relevant to the formation of massive stars in ultradense cores of stellar 
clusters. We find that radiative feedback is more effective than in the case of previous 
calculations in spherical symmetry. The Bondi-Hoyle geometry implies that material 
is flowing nearly tangentially when it experiences the sharp radiative impulse at the 
dust destruction radius, and consequently it is readily perturbed into outflowing or- 
bits. We find that it is difficult for stellar masses to grow beyond around 10 (for 
standard interstellar dust abundances). We discuss the possible implications of this 
result for the formation mechanism of OB stars in cluster cores. We end by proposing 
a series of conditions which must be fulfilled if Bondi-Hoyle accretion is to continue. 
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1 INTRODUCTION 

Massive stars are immensely important astrophysical ob- 
jects, far more so than their rarity would imply. They are 
incredibly bright, being visible at cosmological distances. 
They manufacture most of the metals in the Universe, both 
during their brief lives, and in the supernovae which mark 
their deaths. Massive stars also exert large forces on their 
surroundings, through winds, ionisation and radiation pres- 
sure. 

Despite this importance, massive star formation is a 
poorly understood process. Observational studies are ham- 
pered by the distance to massive star forming regions, and 
the high degree of obscuration in such regions. From a the- 
oretical point of view, the very existence of massive stars 
presents a challenge. Early studies (in spherical symmetry) 
suggested that radiation pressure on dust should prevent 
the assembly of such stars unless the opacity of the dust 
is arbitrarily reduced by an order of magnitude compared 
with interstellar values JWolfire fc CassineUil Il986l. Il987ll 
More recently, simulations under conditions of axisymme- 
try have shown that accretion via a disc c an ease the prob- 
lem somewhat jYorke fc SonnhalteJliool . However, other 
authors have continued to seek other possible formation 
me chanisms for m assive stars. One suggest ion ( postulated 
bv iBonnell et all (il998) and developed as iBonnell fc Batei 
112003)) is that massive stars are formed in the ultra-dense 
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cores of young stellar clusters. Such core s are associated with 
OB star formation (iBlum et alJl20o'3l . alth ough observed 
cores do not reach the densities required bv IBonnell et all 
In the merger models of BonncU ct al., continued accretion 
onto the stars drives an adiabatic contraction of the c luster 
core to the high densities required. IBonnell et al.l i200ll) used 
hydrodynamical simulations of cluster formation to show 
that accretion onto the stars in the core is reasonably well 
described by Bondi-Hoyle accretion (i.e. the accretion cross 
section is appropriate to a compact object moving super- 
sonically through a gaseous medium). However, such sim- 
ulations include no feedback, and thus cannot address the 
question of whether massive stars can continue to accrete by 
this mechanism in an ultra-dense cluster core. 

The interaction of outflowing radiation with a dusty ac- 
cretion flow may be envisaged in the following terms: Dust 
sublimes at a radius in the flow, r^j, where its temperature 
attains a value of 1700 K. Inward of this radius, there is 
no significant momentum exchange between the radiation 
field and the gas, until very small radii, where the gas be- 
come ionised. Radiation that impacts the dusty flow at r^ is 
absorbed within a very narrow region, since the Rosseland 
mean opacity of the dust at the stellar radiation temper- 
ature (~ 30000 K) is very high. The absorption of stellar 
luminosity L in this region is associated with momentum 
transfer to the gas at a rate L/c and the inflow suffers a 
sharp deceleration at this point. This absorbed radiation is 
subsequently thermalised, i.e. re-emitted at a temperature 
equal to the local dust temperature (~ 1700 K). The rate of 
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momentum transfer to the flow that is associated with this 
thermaUsed field is ~ t^ L/c, where tr is the optical depth 
of the flow to the thermalised radiation. 

In this paper, we explore the issue of radiative feedback 
by radiation pressure on dust in the Bondi-Hoyle geome- 
try. In Bondi-Hoyle accretion of isothermal gas, the flow is 
gravitationally focused behind the massive object (star) and 
forms a shock in the downstream direction. For gas enter- 
ing with an impact parameter less than a critical value (the 
Bondi-Hoyle radius), the loss of tangential momentum ren- 
ders the shocked gas bound to the star and enables it to 
be accreted by way of radial infall along the shock. Such ac- 
cretion column geometry evidently has different implications 
for the interaction between the mass flow and radiation field 
than either of the geometries in which radiative feedback has 
been simulated to date (i.e. spherical or disc accretion). Our 
study will be a first look at how radiative feedback on dust 
operates in Bondi-Hoyle geometry We will also be mind- 
ful that in the environment in which such accretion may 
be important, it is not only the geometry that is different 
from conventional studies but also the physical parameters. 
In particular, the hydrodynamical simulations suggest that 
the mean gas densities in the hypothesised ultra-dense core 
phase is extremely high (10® — 10^ cm~'^) and we will find 
that this is an important factor in determining the efficacy 
of radiative feedback. 

In the simulations described below, we employ an 
algorithm for simulating the effect of radiative feedback 
on a dusty flow that we developed in a previous paper 
jEdear fc Clarkd lioo:^ . The aim of this algorithm is to pre- 
serve the sharp deceleration at the point where the dust 
sublimes, while avoiding the co mputational cost of a f ull fre - 
quency dependent calculation. IWolfire &: Cassinelli (^S^) 
showed that a frequency averaged (grey) approach missed 
this sharp impulse (which correspon ds to L/c worth of mo- 
mentum) . The reader should refer to lEdear fc Clark3 ll2003l) 
for a full description (and test) of our method, but the gen- 
eral algorithm is as follows: 

• The radiation field is split into direct and thermalised 
components 

• The direct field is attenuated using wavelength depen- 
dent opacities 

• The balance of the energy is placed into the thermalised 
field, which is then solved using the diffusion approximation 

The initial location of the dust destruction front is estimated 
based on the temperature of a bare grain in the stellar ra- 
diation field, but the algorithm iterates if the thermalised 
field implies the dust melts further from the star. The me- 
chanical effect of anisotropic scattering is included in the 
radiation pressure opacity, but the scattered radiation is 
not followed further. In this paper, we solve the diffusion 
approximation along each radius under the assumption of 
spherical symmetry. This was done for reasons of simplicity, 
although it should be possible to relax this approximation 
in future work. Spherical symmetry is obviously not a good 
approximation for the accretion column, but we shall show 
that the fiow is disrupted in a region where the deviations 
from spherical symmetry are not large. This gives us some 
confidence that our conclusions are not merely an artifact 
of our algorithm. 

The structure of this paper is as follows; In section|5|we 



Figure 1. Sketch of the Bondi-Hoyle accretion geometry 

present a brief mathematical analysis of Bondi-Hoyle flow. 
Section 13 details our search for suitable simulation parame- 
ters. The results from our hydrodynamic simulations follow 
in section 2] We discuss our results and present out conclu- 
sions in sections 13 and |S| respectively. 



BONDI-HOYLE ACCRETION: ANALYTIC 
APPROACH 



iHovle fc Lvttletoiil (Il939ll considered accretion by a star 
(mass M) moving at a steady speed {v^) through an in- 
flnite gas cloud. The gravity of the star focuses the flow into 
a wake which it then accretes. The geometry is sketched in 
flgure 

By considering ballistic orbits, iHovle fc Lvttletoiil con- 
cluded that all material with an impact parameter satisfying 



, , 2GM 
C < Chl = 



(1) 



would be accreted. If the cloud has density p^, the mass 
flux is 



(2) 



which is known as the Hoyle-Lyttleton accretion rate. By 
considering the stability of the accretion co lumn (the wake 
following the point mass on the ^ = axisl . iBcindi fc Hovld 
(1944) concluded that the true accretion rate could be as 
littl e as ha l f this value . 

iBondi 1^53) studied spherically symmet ric accretion 
onto a point mass. The analysis shows (see e.g. iFrank et alJ 
i2002l) ') that a Bondi radius may be defined as 



2GM 



(3) 



where c^b is the sound speed at vb ■ Flow outside this radius 
is subsonic, and the density is almost uniform. Within it, 
the gas becomes supersonic and moves towards a free fall so- 
lution. The similarities between eguationslTlandl^led lBondil 
to propose an interpolation formula: 



M = 



,,2 ^3/2 



(4) 



This is known as the Bondi-Hoyle accretion rate The cor- 
responding is formed by analogy with equation p] On 
the basis of their numerical calculations. IShima et all (ll985^ 
suggest that equation^Jshould acquire an extra factor of two, 
to become 



(c^ 



)S/2 



(5) 
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Figure 2. Sketch of conditions at the dust destruction radius destruction front is marked by the dotted Hne 



which then matches the original Hoyle-Lyttleton rate as the 
sound speed becomes insigni fi cant. 

iBisnovatvi-Kogan et alJ 1^7^ derived the foUowing 
equations to describe the flow, based on the ballistic ap- 
proximation: 



2 2GM Cv'i 
'"L H T- 



ve 



r 



GM( 1 + cos 6*) +C«2o sine 

PooC 



(6) 
(7) 
(8) 
(9) 



The first three equations are obtained from the Newtonian 
orbits. To obtain the final equation, one must solve the gas 
continuity equation using the previous three. 

2.1 Simple Feedback 

The simplest way of including radiative feedback is to place 
L/c worth of momentum into the flow at the dust destruc- 
tion radius, r^. This misses the effect of the thermalised 
radiation field, but is a useful first approximation. 

Consider the flow sketched in figure|21 The 'd' subscripts 
denote quantities immediately outside r^, 'i' subscripts de- 
note values immediately inside r^. All the velocities are ra- 
dial, since the radiative impulse is radial. The velocities 
will be unchanged. The flow is incident from the left, en- 
counters the dust destruction front, and changes density and 
velocity. Assuming that the gas remains isothermal, we can 
conserve mass and momentum to find 



27i"Pd*^d'"d sin ^ d^ — — sin 6 d6 
2c 



(10) 
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where the sin6'd6' terms are taking care of solid angles. 
This pair of equations bears a suspicious resemblance to the 
Rankine-Hugoniot conditions for shocks. We can eliminate 
the value of between these equations, to give 

L 



■ Vd 



(12) 



At first sight, this appears to be independent of angle. How- 
ever, it is not - and are somewhat complex functions 
of 0. 

In simple physical terms, the flow will slow dramatically 
in the radial direction as it passes through r^. This will cause 



gas streamlines to move onto more circular orbits (since their 
angular momentum is unaffected). The gas streams could 
then flow around inside r^, and still end up on the accretion 
column trailing the star. Indeed, if were larger than Cbhi 
then a small radiative impulse could increase the accretion 
rate. Material originating from large ( would be moved onto 
orbits closer to circular (since the radial impulse would not 
affect angular momentum), and hence be more likely to be 
bound to the star after encountering the accretion column. 
For a protostar in a protocluster, this would require > 
36kms"^, which is a little on the high side. The accretion 
rate (proportional to - cf equation ^ would also be 
rather low, unless the density were increased to compensate. 

There will be a gas streamline which has its closest ap- 
proach to the star at r^. Eauation ll2l is obviously not going 
to apply to this streamline, since it will have = 0. If 
the angle, 6c, at which this occurs is close enough to zero, 
this will not be a problem, since it will be lost in the wake. 
However, there will be a range of angles with 9 < 9^ for 
which equation 1121 will imply a negative value for u;. Phys- 
ically, this corresponds to the flow being turned around at 
r^. While this is fine for particles, it will be disastrous for 
the ballistic approximation of the gas. Streamlines will cross, 
destroying the assumptions necessary to treat the gas using 
simple Newtonian theory. 

Using equation 1121 and the fact that the velocity in the 
6 direction is unchanged, it is fairly straightforward to com- 
pute the trajectory followed by the gas inside the r^ (on the 
assumption of ballistic flow) . 

Figure |21 plots some sample trajectories for flow past 
a 10 M© star. In this plot, = 10~^® gcm~"^ and = 
5 X 10^ cms"^. We shall show in the next section that these 
parameters are appropriate for the ultra-dense cores postu- 
lated by BormcU & Bate. To make the perturbations to the 
flow visible, we used an artificially large value for the lu- 
minosity of 10** Lq (approximately twice the ZAMS value 
for a 10 Mq star). We placed the dust destruction radius 
at r^j = 5 X lO^** cm, which is a reasonable value for these 
conditions. 

Our simplified feedback model implies that radiation 
pressure on dust will give rise to orbit crossing (and hence 
the formation of a shock cone) in the downstream direction 
(i.e. for 6 < 9^). For a luminosity appropriate to a 10 M© 
ZAMS star, 9^ ~ 0.27r rad, which far exceeds the width of 
the unperturbed Bondi-Hoyle wake. 

To proceed further, it is necessary to resort to full hy- 
drodynamic calculations, both because these can follow the 
flow in the orbit crossing regime and because they can in- 
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elude additional feedback from the thermalised radiation 
field. This important effect, omitted in the above analysis, 
decelerates the flow before it encounters the dust destruc- 
tion front. In addition, we can use hydrodynamical simu- 
lations to explore the case that the equation of st ate of 
the gas is n ot isothermal, since it is well known (e. g. iH untl 
Jl97ll I1979II and Ruffcrt ( 19941: iRuffert fc ArnetTi 199411 ) 
that Bondi-Hoyle flow is no t well described by the equations 
of .Bisnovatvi-Kogan et alJ in this case. 



3 PARAMETERS 

Simulating Bondi-Hoyle flow in the context of star forma- 
tion is rather difficult, due to the large range in length 
scales which must be resolved. Although we are interested 
in stars rather more massive than the compact objects typ- 
ically studied in Bondi-Hoyle simulations, the protostellar 
velocities in a cloud core are generally much lower as well - 
typically only a few kms"'^. This tends to make ^hl rather 
large (cf equation 0, and many times larger than the proto- 
stellar radius. We can reduce the problem, since we are not 
interested in following the flow right down to the protostar, 
but only to the dust destruction front, r^.^ However, the 
ratio Chl/'^cI can still be painfully large. Another motivation 
for keeping (^jjl small is self-gravity. If the mass enclosed in 
the simulation is comparable to that of the protostar, consis- 
tency requires that the flow is self-gravitating - which would 
require a great deal of computational effort. The Hll region 
surrounding the protostar must also be kept away from the 
dust. If ionised gas reaches the dust, then it might be able 
to destroy the dust before the grains reach their nominal 
sublimation temperature. 

To find suitable simulation parameters (M, v^, poo), we 
make use of the analy tic formulae for pressure f ree B ondi- 
Hoyle fiow derived by iBisnovatvi-Koean et alJ lll979ll . We 
estimate the volume of ionised gas by balancing the num- 
ber of ionisations and recombinations along each sight line, 
in a straightforwa rd generalisation of the Stromgren Sphere 
JStromgrenlll939ll . To calculate the number of io nising pho- 
tons, we make use of the ZAMS formulae of iTout et alJ 
il996l) and the Stefan-Boltzmann Law. We seek parameters 
which satisfy the following conditions: 

• The predicted accretion rate will cause a significant 
change in mass on a timescale of a million years or so 

• The volume of ionised gas must lie within 

• The ratio r^/CBH rnust not be too small 

Note that the first two requirements are motivated by physi- 
cal considerations. The third is a computational practicality. 
Unfortunately, these requirements tend to conflict with each 
other - a faster gas flow will help reduce CbH) but it will also 
reduce the central concentration of the fluid. This will let 
the ionised region grow. 

Base d on information supp lied by Ian Bonnell (pub- 
lished as lBonnell fc Bat^ (|200j) the reasonable ranges for 
the relevant parameters are given in Tabled However, the 
upper end of each range is beginning to push limits to the 

^ is calculated by computing the sublimation temperatu re of 
a dus t grain in the stellar radiation field. See .Edgar & Clarkd 
^2003^ for details 



Gas Density 10~^* gcm""^ < < 10~^^ gcm~^ 
Gas Velocity 10* cms"^ < < 10^ cms"^ 

Stellar Mass 10 Mq < M < 100 Mq 

Table 1. Range o f Parameters for an ultra-dense core 
iBonnell fc Batdl2002l) 

Gas Density 10"^'' gem"'' 
Gas Velocity 5 X 10^ cms"^ 
Stellar Mass 10 Mq 

Table 2. Chosen parameters for Bondi-Hoyle simulations 

extreme. In particular, there is little point seeking Bondi- 
Hoyle accretion solutions for a 100 Mq star, unless we have 
already established that it can grow to this mass. Having 
said this, more massive stars make it easier to keep the ioni- 
sation problem to a minimum. This is because is an 'area' 
effect (cf the Stefan-Boltzmann Law) whereas ionisation is 
a 'volume' effect (cf the Stromgren Sphere). 

From this, the parameters given in Table |21 squeeze 
through the selection criteria. These give an accretion rate 
of 3 X 10"* M0 yr"^ and (bh = 1-1 x 10^*^ cm. The dust de- 
struction front due to direct stellar radiation is located at 
r^j = 2.1 X lO'^* cm. The ionised gas should always lie within 
~ 1.3 X lO'^* cm of the star,^ which is sufficiently inside 
for our purposes. This is helped by the following considera- 
tions: 

• The value of is likely to increase, due to the ther- 
malised radiation field 

• The escape velocity at the edge of the Hll region is on 
the order of forty kilometres per second, while the sound 
speed in the ionised gas is only a few kilometres per second 

Therefore, ionisation is unlikely to be a problem at the start 
of the simulations. So long as accretion continues, the Hll 
region is unlikely to grow. However, if the radiative feedback 
causes significant disruption to the fiow, then falling densi- 
ties within are likely to allow the ionised volume to grow, 
violating our assumptions. 



4 DYNAMICAL SIMULATIONS 

The Bondi-Hoyle solution does not occur for gaseous fiow. 
Shocks and instabilities occur, requiring numerical solutions. 
To simulate the ffow, we used the lateste d publically avail- 
able v ersion (2.0.3) of the Zeus2d code of lStone fc NormanI 
il992l) in spherical polar mode. Radiative feedback was in- 
cluded using the radia tive transfer algorithm developed by 
lEdear fc Clark3 ll2003h . with the radiative feedback appear- 
ing as an extra source of momentum. In the present work, 
we solved the diffusion approximation (and hence as well) 
along each ray using the assumption of spherical symmetry. 
This should be a good approximation on the upstream side, 
but will be poor in the wake, where the density gradient will 
tend to transport radiation sideways. This would have the 
effect of expanding the wake, and reducing its ram pressure. 
Despite the computational efficiency of our new algorithm, 

^ Note that this is still subst antially small er than the Hypercom- 
pact Hll regions discussed by iKurta l2002l) 
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it is still expensive to solve on every timestep. Accordingly, 
the density structure along each ray was stored whenever the 
radiative feedback was calculated. The solution for a partic- 
ular ray was only recomputed when this structure changed 
signficantly. 

We calculated accretion by depleting the innermost grid 
cells of material, by a fixed fraction each timestep. On the 
upstream side, the boundary conditions were those given by 
equations|S]to|3] This enabled th e outer boundary o f the sim- 
ulation to be moved inwards fcf lKoide et al.l lll99ll) ). easing 
the computational load. On the downstream side, we im- 
posed outflow conditions. For the equation of state, we tried 
7 = 1,4/3,5/3. Although the d ensities suggest isothermal - 
ity is the best approximation JLow fc Lvnden-Belll[l973) . 
there will also be radiative heating and possibly other pro- 
ce sses. T he star was parameterised by the ZAMS formulae 
of lTout e t al. (1996). 

In the majority of our calculations, we used 50 grid 
cells in the 8 direction, and 160 in the radial direction. The 
angular spacing was uniform, but the radial spacing was 
not. Only the inner five radial zones (extending to r = 0) 
were uniformly spaced, and formed the accretion cavity. The 
remaining 155 grid cells were logarithmically spaced, joining 
smoothly on to the uniform grid. For some test cases, we 
increased the resolution, but this did not change our results 
significantly. The resolution we used means that most of the 
direct stellar impulse is applied in a single grid cell, but 
Zeus2d does not seem to mind this (cf section [4.21 . In the 
isothermal case, the accretion column is covered by about 
five azimuthal zones. 

In all our calculations, the only source of gravity is the 
star itself. The self-gravity of the gas is neglected, due to 
the computational cost. However, the neglect is not with- 
out justification. The gas mass contained in our simulations 
should remain less than the stellar mass. Furthermore, these 
simulations are meant to be a small sample around a larger 
cluster core. In this case, the structure of the cluster itself 
would provide some support against local gravitational col- 
lapses of the gas. As we shall see, isothermal gas can give 
rise to high gas densities in the accretion column, and it is 
possible that self-gravitating clumps could form. However, 
such clumps would still be accreted by the star. Finally, any 
over-densities which do form will tend to absorb more radi- 
ation, and the increased temperature will tend to stabilise 
them against gravitational collapse. 
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Figure 4. Accretion rates for plain Bondi— Hoyle flow 
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Figure 5. Density contours for high Mach number flow around 
a smaU accretor 



4.1 Flow without Feedback 

As mentioned above, the ratio CuL/'^d is rather large for our 
simulations - larger than most simulations of Bondi-Hoyle 
flow to date. Therefore, we first simulated Bondi-Hoyle flow 
without feedback, to provide a benchmark for the runs with 
feedback included. 

Figure 2] shows the resultant accretion rates as a func- 
tion of time. The simulations with 7 values of 4/3 and 5/3 
settle into a quasi-steady state, with the stiffer equation of 
state giving a slightly lower accretion rate. The isothermal 
run shows quite strong oscillations in its accretion rate, al- 
though the mean value is close to Mhl- These oscillations 
are likely to be due to the gas suddenly shocking as it en- 
counters the ^ = axis. In a very short distance, it has to 
turn around (the test of accretion in the standard analysis 



is for material flowing away from the accretor), and it is 
evidently unable to do this in a steady manner. 

We plot sample density and velocity fields for these 
runs in Figures |5] and |S] Note that the size of these plots 
is rather smaller than i^hl ~ 10^^ cm. While the isothermal 
ru n is obviously quali tatively similar to the flow described 
bv lHovle fc Lvttletonl the stiffer equations of state lead to 
rather different behaviour. For 7 = 5/3, a bow shock forms, 
with the downstream flow close to spherical (cf e.g. iHund 
iigTll) *). Putting 7 = 4/3 moves the bow shock back, attach- 
ing it to the accretor. The isothermal flow is qualitatively 
similar to that described by equations|S|to|^ We found that 
the oscillating accretion rate plotted in Figure 01 was due to 
a slight expansion and contraction in the accretion column, 
rather than a dramatic change in the flow pattern. 
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Figure 6. Velocity fields for high Mach number flow around a 
small accretor. The approximate location of the shock front is 
shown by the dotted line. Flow is incident from the left 



We can also make a useful comparison between the ram 
pressure in the radial direction and the sharp impulse at the 
dust destruction radius. At the most basic level, accretion 
can only proceed if 



> 1 



(13) 



(the angular factors are the same for both ram and radia- 
tion p ressure). For a 10 M© star, the formulae of lTout et alJ 
l)l99fil predict L = 5552 . Results for an expected dust 
destruction radius of 2.1 x 10^'' cm are plotted in Figure |7| 
From this, it is obvious that we have happened upon an in- 
teresting region of parameter space. The adiabatic runs have 
H values which get as low as a few. This means that the di- 
rect impulse will not be sufficient to halt the flow. However, 
the extra momentum from the thermalised field may slow 
the flow sufficiently to allow the accretion flow to be shut 
down. The S values for the isothermal case vary even more 
dramatically. Just outside the wake, the radiation pressure 
at the dust destruction front will dominate the ram pressure 
(unsurprising, given the velocity held shown in Figure 0. In 
the wake, the ram pressure is orders of magnitude greater 
than the radiation pressure. 



4.2 Direct Impulse Only 

To provide comparison with later simulations, we per- 
formed three runs with only the direct impulse included. 
We achieved this by injecting L/c worth of momentum into 
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Figure 7. Ratio of radial ram pressure to L/c for high Mach 
number flow around a small accretor 
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Figure 8. Accretion rates for runs with the direct impulse only. 
The corresponding accretion rates from figure |1] arc plotted for 
comparison 



the grid at a fixed radius. As may be seen from figure |S1 
the effect on the accretion rate is minimal. The greatest dif- 
ferences occur for the isothermal case, and even then, the 
changes probably have more to do with the intrinsic insta- 
bility of the flow, rather than any significant changes due to 
feedback. These results are not too surprising in the light 
of figure which showed that the ram pressure should ex- 
ceed the radiation pressure for most of the flow at the dust 
destruction radius. 
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Time / Crossing Times 

Figure 9. Accretion rates for runs with reduced dust. The corre- 
sponding accretion rates from figure |3 are plotted for comparison 

4.3 Feedback with Reduced Dust 

Three ru ns were performed w i th the depleted dust abun- 
dances of IWolfire fc Cassinellil il987ll . This involves an ad 
hoc reduction of the dust abundances b y a factor of ten, and 
also a reduction in graphite grain size. IWolfire fc Cassinelli 
found this reduction was necessary to ensure material was 
always accelerating inwards in their calculations. 

Accretion rates for these runs are plotted against time 
in figure |5| From these, it is apparent that the addition of 
feedback is not having a dramatic effect on the flow. The 
run with 7 = 5/3 now shows larger fluctuations in its ac- 
cretion rate, and the overall rate has dropped too. Smaller 
changes have occurred for the 7 = 4/3 gas, with the accre- 
tion rate only falling slightly when compared to the pure 
hydrodynamic run. For the isothermal case, the change in 
accretion rate is barely visible compared to the large fluctu- 
ations present. 

On examining the density contours, we found that the 
run with 7 = 4/3 was little changed - unsurprisingly, in 
view of Figure |3] The isothermal simulation showed a small 
density enhancement adjacent to the wake, corresponding to 
the portion of Figure|7|where H was less than unity. Beyond 
this, there were few changes. There were greater changes 
for the 7 = 5/3 gas - the flow was far less smooth (as one 
might expect, given the changing accretion rate). Since the 
reduced dust abundances made the Rosseland optical depth 
low, the shape of the dust destruction front was close to 
spherical in all cases. The only signiflcant deviation was for 
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Figure 10. Comparison between the Rosseland optical depth (t) 
and the ratio of ram and radiation pressures (H) for the isothermal 
run with full dust 

the accretion column in the isothermal case. Even here, the 
dust destruction front was only pushed out slightly. 

4.4 Feedback with Full Dust 

Increasing the dust abundances to their normal galactic val- 
ues had a dramatic effect on the flow. Soon after the impulse 
switched on, material ceased penetrating r^, and the grid 
cells inside drained of material. Zeus2d was then forced to 
terminate. Consider Figure [Till This plots H and tr values 
for the isothermal gas shortly after the impulse is switched 
on. If the gas is outflowing, H is set to zero. It is relevant 
to compare these two quantities because the amount of mo- 
mentum imparted to the flow by the thermalised field is 
Ltr/c. Since H measures the ram pressure of the flow at 
in units of L/c, if tr ^ H, then the thermalised radiation 
should reduce the ram pressure at enough to allow the 
direct radiation field to halt inflow. Figure 1101 shows that 
the total momentum imparted to the flow by the luminosity 
completely dominating the ram pressure in the flow. This is 
the reason accretion halted. 

When we examined the position of the dust destruction 
front, we found that it was still close to spherical for most of 
the flow. However, the attached shock in the 7 = 4/3 case, 
and the accretion column in the isothermal case caused no- 
ticeable deviations. In the isothermal run, the dust destruc- 
tion front in the accretion column was located roughly four 
times further from the accretor than on the upstream side. 

We attempted various numerical 'fixes' in order to en- 
sure that numerical problems were not interfering with the 
physical analysis. In some runs, we spread the effect of the 
direct stellar radiation field over a number (typically seven) 
of grid cells. In others, we tried adjusting the CFL condition, 
to incorporate a term based on s = O.bat^ . Finally, we tried 
both methods together. None of these 'worked' - accretion 
was still halted in every case. We therefore concluded that 
our results were not an artefact of the numerical method. 
This is consistent with the theoretical arguments discussed 
above. 



5 DISCUSSION 

Radiative feedback has a much more dramatic effect on 
Bondi-Hoyle accretion than spherically symmetric accre- 
tion. This is brought about by the lower central concen- 
tration and presence of centrifugal support found in Bondi- 
Hoyle flow. 



8 R.G. Edgar & C.J. Clarke 



Although we have only been able to study one set of 
parameters, there are several implications for the formation 
of massive stars in clusters. Since L ~ for stars close to 
10 M0 , the fact that our simulations were obviously on the 
borderline between accreting and not accreting implies that 

• Less massive stars will have few problems accreting in 
the Bondi-Hoyle geometry 

• More massive stars are unlikely to accrete significant 
amount of mass under these conditions 

Unfortunately, it is not easy to determine an exact limit - 
especially since we have not been able to test a wide range 
of parameters, and the fact that it is the reprocessed radi- 
ation which is of critical importance. Consider eauation ll2l 
Since L/r^ should be approximately constant (cf the Stefan- 
Boltzmann equation), if the 'inside' radial velocity is posi- 
tive (indicating accretion) for a protostar at some mass, it 
should remain positive as the mass increases. The sharp im- 
pulse at the dust destruction radius cannot stop the inflow 
by itself. It can only do so if the flow has been slowed by 
the thermalised radiation. Of course, if the value of rises 
too much (due to the contraction of the cluster), then the 
accretion timescale (cf equation |2l will become longer than 
the main sequence timescale of the star. 

Estimating the effect of the thermalised radiation is 
much more difficult, since the impulse from it is distributed 
throughout the flow. Based on the findings presented above, 
requiring that H > tr,"^ at the dust destruction radius is the 
most reasonable condition for the purpose of making esti- 
mates. However, this should be treated with caution, since 
it is not rigorously derived. 

Despite these problems, we can still suggest regions of 
parameter space which might permit accretion, and hence 
be worthy of future investigation. These regions must fulfil 
several conditions: 

• -^HL must be large enough to permit the star to gain 
appreciable mass before it leaves the main sequence 

• They must have H > 1 at the dust destruction radius 
(or the flow will be reversed) 

• The Hll region must lie within the dust destruction ra- 
dius 

• They must have E > tr at the dust destruction radius 
(again, to avoid reversal of the flow) 

For a given mass of star, these conditions define a region of 
(Poo,t^oo) space which are likely to permit accretion in the 
Bondi-Hoyle geometry. In Figure 1111 we show the allowed 
parameters for a variety of stellar masses, evaluated using 
equations |S| to |n| for the full dust abundances. For the first 
condition, we required A/jjl > 10~^ M0 yr. The last three 
conditions above are functions of angle, so in these cases we 
tested the condition on the upstream [6 — tt) direction. The 
Hll region will be largest in this direction. Furthermore, the 
analytic approximation will hold to smaller radii. We also 
added the condition that ^hl > Sr^, and took the edge of the 
calculation to lie at 2^jjl (for the purposes of calculating tr). 
Both of these conditions are somewhat arbitrary. First, note 
that the parameters tested above {p^ — 10^^® gcm~'^ and 

^ Recall that H is the ratio of ram to radiation pressure at - 
see equation 1131 



10' 



105 


















- 1 1 


Ms 


















• ^ 




- 10 


Mq 










-■ 9 


Me 








: 


-- 8 


Mg 






: 




— 7 


Mo 






















yy' -■ 














: / 

/' 










yy y 












■ ^ 


/ 












T" , .L— T^-r^rT , 1 






-18 




10-" 






10 



p„ / g cm 3 

Figure 11. Allowed parameters for Bondi-Hoyle accretion for a 
variety of stellar masses, as suggested by the conditions in the 
text. The contours enclose the allowed regions for the marked 
stellar masses, when the dust abundance is assumed to be Galac- 
tic. The star marks the location of the simulations presented in 
this paper 

Woo = 5 X 10^ cms'^ for M = 10 M©) lie outside the allowed 
range for the full dust (this point lies inside the permitted 
region, as calculated for the reduced dust). Of course, this 
is not a full test of the validity of these conditions, but it 
is encouraging. The sloping straight line boundary is caused 
by the requirement on M . The horizontal boundary which 
appears close to = 2 x 10^ cms"^ for the higher masses 
is set by the constraint on ("hl- 

For accretion to occur. Figure [TTI suggests that the op- 
tical depth must be kept low, either by high velocities (so 
that the distance to the outer boundary is small) or low 
densities. The permitted region for the higher mass stars at 
high velocities is almost independent of density because p 
does not vary much between and 2("hl- This means that 
p cancels from the ratio H/tr, leaving the ratio a function of 
Woo only. At low densities, the rather complicated shape of 
the allowed region is due to the solution of the thermalised 
radiation field, and the calculation of tr. Note that there is 
no allowed low density region for the 11 Mq case. 

It seems that Bondi-Hoyle accretion is going to run 
into serious trouble at about 10 M©. Even for masses of 
only 8 Mq , the allowed region is curtailed quite severely. 
However, these are only estimates, and assume that Bondi- 
Hoyle accretion will be a good approximation over all the 
parameter space covered by figure ITTl Unfortunately, this is 
unlikely to be the case. For example, a velocity of lO** cms"^ 
for a 10 Mq star gives i^hl '~ 10^^ cm, which is far larger 
than the interstellar separation in a cluster core.* Also, 
a 10 Mq star moving at 10"" cms"^ through at cloud with 
Poo ~ 10~^^ gem""' would have an unreasonably high ac- 
cretion rate of Mhl ~ O.SMQyr"'^. These concerns would 
restrict the allowed parameter space further. 

There are further problems, which our simulations do 

^ Additionally, simulating such a large volume would be compu- 
tationally painful, due to the enormous difference in length scales. 
We will still have ~ 10^* cm. Worse still, a cloud that big is 
going to be self-gravitating 
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not address. We have not included an accretion luminosity, 
since it will complicate the feedback mechanism. However, 
the accretion luminosity implied by the parameters of Ta- 
ble |5| is roughly three times the ZAMS luminosity. A short 
test, on the assumption that the accretion luminosity was 
also isotropic (this may not be the case), suggested that 
accretion was not possible even for the reduced dust abun- 
dances. Our inability to simulate the ionised region is also a 
concern. Flow in the Bondi-Hoyle geometry is far less cen- 
trally concentrated than in the spherically symmetric case, 
so the Hll region might be able to grow. If it can reach 
the dust destruction radius, some rather unpleasant chem- 
istry would have to be simulated as well. Our imposition 
of axisymmetry and low resolution in the inner portions of 
the grid also prevent us examining the question of accretion 
discs. If the flow at infinity is non-uniform (rather likely in a 
protocluster) , then there is the potenti al for angular momen- 
tum accretion too iRufrertll99lll999t) . This would drive the 
formation of an accretion disc. 

The problems radiative feedback presents in the Bondi- 
Hoyle geometry have potentially serious implications for 
the merger model o f massive star formation proposed by 
iBonnell et al.l il998tl . Accretion drove the contraction of 
thei r clusters to the high st ellar densities required for merg- 
ers. iBonnell fc Batd l|2003) presented a full hydrodynamic 
simulation of such a cluster, but did not have a cut off mass 
for accretion (their model was scale free). If the most mas- 
sive stars in the cluster core stopped accreting, then it is 
possible that the cluster core would not contract enough to 
drive a high merger rate. However, this is very uncertain. 
IBonnell et al.l ||2001j) simulated another cluster, and their 
results show many low mass stars in the core, in addition 
to higher mass ones. It is possible that the lower mass stars 
could accrete enough material to drive the stellar density 
to tha t required for mergers. Note also that IBonnell et alJ 
il998h did prevent stars more massive than 10 accreting 
in their simulations, but still found merging stars. 

6 CONCLUSION 

We have shown that if stars form in the ultra-dense core s of 
young clusters, as postu l ated i n the models of^onnoll ct alJ 
lil99d) :lBonncll fc Bat j J2002l) . then the feedback from ra- 
diation pressure on dust is likely to prevent accretion on to 
stars more massive than ~ 10 Mq . This feedback is more 
disruptive to the accretion flow than in our previous sim- 
ulations (Edgar & Clarke 2003) which used the initial con- 
dition s employed in the classic work of IWolfire fc Cassinelli 
l|l987fl (i.e. spherically symmetric accretion and initial gas 
densities of 10"^^ gcm"^). 

The reasons for this more effective feedback in the 
present case are two-fold. Firstly, the ul tra-dense cluster 
cores that are invoked by IBonnell et alJ are unobservably 
short lived. In this phase the local gas density of their cores 
is extremely high, exceeding that s een in regions of mass ive 
star formation, or that employed bv lWolfire fc CassinellA by 
many orders of magnitude. Consequently, the optical depth 
to the thermalised radiation field (tr) is very high. Since 
the total rate of momentum input into the gas is ^ tkL/c, 
this increases the coupling between the radiation and gas 
dramatically. Secondly, in the case of spherically symmetric 



accretion, the gas density is centrally concentrated and the 
velocity field is radial. Both of these effects enhancing the 
radial component of the ram pressure in the accretion flow 
at the dust destruction radius. In the case of cluster cores, 
by contrast, accretion proceeds in Bondi-Hoyle geometry, 
and flow is nearly tangential as it approaches the accretion 
wake (see Figure 0. Consequently the radial component of 
the ram pressure is relatively small, and this again favours 
the disruption of the accretion flow. 

The above results imply that if massive stars indeed 
form in these hypothetical ultradense cores, then for masses 
greater than ~ 10 M0 , they must be assembled entirely by 
stellar collisions and consequent mergers rather than by ac- 
cretion. The attainment of these ultradense conditions is 
however driven by accretion by stars in the core and it cur- 
rently unclear whether accretion on to low mass stars (i.e. 
those with mass < 10 M0) will suffice for this purpose. This 
issue needs to be explored by further simulations that incor- 
porate feedback from radiation pressure on dust into realistic 
cluster simulations. 
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